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Abstract. Given a homotopy connecting two polynomial sys- 
tems we provide a rigorous algorithm for tracking a regular ho- 
motopy path connecting an approximate zero of the start system 
to an approximate zero of the target system. Our method uses 
recent results on the complexity of homotopy continuation rooted 
in the alpha theory of Smale. Experimental results obtained with 
the implementation in the numerical algebraic geometry package 
of Macaulay2 demonstrate the practicality of the algorithm. In 
particular, we confirm the theoretical results for random linear ho- 
motopies and illustrate the plausibility of a conjecture by Shub and 
Smale on a good initial pair. 

The numerical homotopy continuation methods are the backbone of 
the area of numerical algebraic geometry; while this area has a rigor- 
ous theoretic base, its existing software relies on heuristics to perform 
homotopy tracking. This paper has two main goals: 

• On one hand, we intend to overview some recent developments 
in the analysis of complexity of polynomial homotopy contin- 
uation methods with the view towards a practical implemen- 
tation. In the last years, there has been much progress in the 
understanding of this problem. We hereby summarize the main 
results obtained, writting them in a unified and accesible way. 

• On the other hand, we present for the first time an implemen- 
tation of a certified homotopy method which does not rely on 
heuristic considerations. Experiments with this algorithm are 
also presented, providing for the first time a tool to study deep 
conjectures on the complexity of homotopy methods (as Shub 
& Smale's conjecture discussed below) and illustrating known 
-yet somehow surprising - features about these methods, as 
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equiprobability of the output in the case of random linear ho- 
motopy and the average polynomial or quasi-polynomial time 
of the algorithms studied by several authors. 

Our project constructs a certified homotopy tracking algorithm and 
delivers the first practical implementation of a rigorous path-following 
procedure. In particular, the case of a linear homotopy is addressed in 
full detail in Algorithm [T] of Section 3.3 



The structure of the paper is as follows. In Section[T]we introduce the 
basic notations and the general problems that we address. In Section 
[2] we recall the definition of approximate zero, condition number, and 
Newton's method, and equip the space of polynomial systems with a 
Hermitian product. In Section [3] the main problem solved in this paper 
is formulated; we describe a certified algorithm to follow a homotopy 
path. An overview of approaches to finding all the roots of a system is 
presented in Section |4j In Section [5] we give an algorithm to construct 
a random linear homotopy with good average complexity. In Section [6] 
we discuss the implementation of our algorithm. Section [7] demon- 
strates the practicality of computation with the developed algorithm 
and discusses experimental data that could be used to obtain intuition, 
in particular, with regards to a longstanding conjecture of Shub and 
Smale. 
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partially done while the authors were attending a workshop on the 
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program hosted by the Fields Institute, Toronto. We thank that insti- 
tution for their kind support. 



1. Description of the problem 

Let n > 1. For a positive integer I > 1, let Vi — Q[A 1; . . . ,X n ] be 
the vector space of all polynomials of degree at most / with complex 
coefficients and unknowns X\, . . . , X n . Then, for a list of degrees (d) = 
(g?i, . . . , d n ) let V(d) = Vd x x ■ ■ ■ x Vd„- Note that elements in V{d) 
are n-tuples / = (fx, . . . , /„) where fi is a polynomial of degree di. An 
element / G Vu) will be seen both as a vector in some high-dimensional 
vector space and as a system of n equations with n unknowns. 

Problem 1. Assuming f G Vu) has finitely many zeros, find approxi- 
mately one, several, or all zeros of f in C n . 
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It is helpful to consider the homogeneous version of this problem: 
For a positive integer I > 1, let Hi be the vector space of all ho- 
mogeneous polynomials of degree I with complex coefficients and un- 
knowns X , . . . , X n . Then, for a list of degrees (d) = (di, . . . , d n ) let 
"H(d) = "H^i x ••• x %dn- Note that elements in %(d) are n-tuples 
ft, = (hi, . . . , h n ) where hi is a homogeneous polynomial of degree di. 
An element h G H(d) will be seen both as a vector in some high- 
dimensional vector space, and as a system of n homogeneous equations 
with n + 1 unknowns. Note that if C G C n+1 is a zero of h G "H(<i), then 
so is AC, A G C. Hence, it makes sense to consider zeros of h G H(d) as 
projective points C G P(C n+1 ). Abusing the notation, we will denote 
both a point in P(C n+1 ) and a representative of the point in C n+1 with 
the same symbol. Moreover, if necessary, it is implied that the norm 
of this representative equals 1. 

Problem 2. Assuming h G ^H(d) has finitely many zeros, find approxi- 
mately one, several or all zeros of h in P(C n+1 ). 

Let d = maxjc?!, . . . , d n } and T> — d\ ■ ■ ■ d n . Note that d is a small 
quantity, but in general V is an exponential quantity. We denote by 
N+l the complex dimension of H^) and V(d) as vector spaces. Namely, 

There is a correspondence between problems [T] and |2j Given / = 

(fl, ■ ■ ■ , fn) € V( d ), 

f — \ ^ /7* IT ™ 

— / j a a 1 ,...,a n jX l > 

ai + ...+a n <d 4 

we can consider its homogeneous counterpart ft = (ft-i, . . . , ft n ) G 
where 

/, _ \ A „i ydi-{ai-\ ha„) yai v a n 

Hi— a ai,...,a„ A A l • ,,A n ' 

ai+...+a n <d; 

If (771, ... , 7/ n ) is a zero of /, then (1, 771, ... , r/ n ) is a zero ft.. Conversely, 
if (Co, • • • , Cn) e P(C n+1 ) is a zero of ft and Co 7^ then . . . , is 
a zero of /. 

2. Preliminaries 

2.1. Approximate zeros and Newton's method. In general, it is 
hard to describe zeros of / G V(d) or ft G Ti(d) exactly. One may ask for 
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points which are "e-close" to some zero, but this is not a very stable 
concept. The concept of an approximate zero of [21] fixes that gap. 
Given / G V(d), consider the Newton operator associated to /, 

N{f)(x)=x-Df{x)- 1 f{x), 

where Df(x) is the n x n derivative matrix of / at x G C n , also often 
called the Jacobian (matrix). Note that N(f)(x) is defined as far as 
Df(x) is an invertible matrix. We will denote 

i 

N(f) l (x) = N(f)1^N(f)(x) 
namely, the result of I iterations of Newton's method starting at x. 

Definition 1. We say that x G C™ is an approximate zero of f G Vu) 
with associated zero rj G C n if N(f) l (x) is defined for all I > and 

\\N(f) l (x)- V \\<^^, />0. 

The homogeneous version of Newton's method [20] is defined as fol- 
lows. Let h G H id) and z G P(C n+1 ). Then, 

N P (h)(z) = z - (Dh(z) \ z x)- x h(z), 

where Dh(z) is the n x (n + 1) Jacobian matrix of h at z G P(C n+1 ), 
and 

Dh(z) \ z s. 

is the restriction of the linear operator defined by Dh(z) : C" +1 — > C™ 
to the orthogonal complement z x of z. Hence, (Dh(z) l^) -1 is a linear 
operator from C n to z x , and N^(h) (z) is defined as far as this operator is 
invertible. The reader may check that Np(h)(Xz) = XNp(h)(z), namely 
Np(h) is a well-defined projective operator. Note that Np(h) may be 
written in a matrix form 

which is more comfortable for computations. As before, we denote by 
Np(h) l (z) the result of I consecutive applications of Np(h) with the 
initial point z. 

Definition 2. We say that z G P(C ra+1 ) is an approximate zero of 
h G *H(d) w tth associated zero ( G P(C n+1 ) if Nf(h) l (z) is defined for 
all I > and 

d R {N w {h)\z)X)< < ^0 } -, l>0, 
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Here d,R is the Riemann distance in P(C n+1 ), namely 

\(z z')\ 

d R (z,z') = arccos^|0| G [0,7r/2], 

where (•, •) and || • || are the usual Hermitian product and norm in 
C n+1 . Note that <1r(z, z') = cLr(Xz, X'z') for A, A' G C, namely cIr is well 
defined in P(C n+1 ) x P(C n+1 ). The reader familiar with Riemannian 
geometry may check that cIr(z, z') is the length of the shortest C 1 curve 
with extremes z, z' G P(C™ +1 ), when P(C n+1 ) is endowed with the usual 
Hermitian structure (see [21 Page 226].) 

Let / G V(d) and let h G H(d) be the homogeneous counterpart 
of /. In contrast with the case of exact zeros, it may happen that 

z = (zo, ■ ■ ■ , z n ) is an approximate zero of h but still . . . , ^\ is not 

an approximate zero of /. In Proposition [3] we explain how to fix that 
gap. 

2.2. The Bombieri-Weyl Hermitian product. In studying Prob- 
lems [I] and [2j it is very helpful to introduce some geometric and met- 
ric properties in the vector spaces V(d) and H(d)- We recall now the 
unitarily-invariant Hermitian product in 7-L(d), sometimes called Kost- 
lan Hermitian product ([2]) or Bombieri-Weyl Hermitian product ([8]). 
Given two polynomials v,w G Hi, 

V - a ao,-,a„XQ° ■ ■ ■ X" n , 

ao+.-+a n =l 



ao+...+a n =l 

we consider their (Bombieri-Weyl) product 

(v,w) — / I / \ ) a a ,...,a„b ao „„. 

ao+ai+...+a n =l xv 

where ~ is the complex conjugation and 

I \ l\ 



Xa ,...,a n )J a \ ■■•«„! 

is the multinomial coefficient. 

Then, given two elements h = (hi, ... , h n ) and hi = (h[, . . . , h' n ) of 
T-i(d), we define 

(h, ti) = (h, ti 1 )+---+ (hn, h' n ), \\h\\ = (h, h) 1 ' 2 . 

This Hermitian product defines a real inner product in H^d) as usual, 

(h,h') R = Re((h,ti)). 
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We also define a Hermitian product and the associated norm in V(d) 
as follows: Given /, /' G V( d ), let h, hi G H(d) be the homogeneous 
counterparts of /, /'. Then, define 

(f,f) = {h,ti), ||/|| = |N- 

From now on, we will denote by § the unit sphere in %{d) for this norm, 
namely 

S = {h G H {d) : \\h\\ = 1}. 

Note that for solving Problem [2j we may restrict our input systems 
h G Hid) to h G S, for zeros of a system of equations do not change if 
the system is multiplied by a non-zero scalar number. 

2.3. The condition number. The condition number at (h, z) G Hid) x 
P(C n+1 ) is defined as follows 

MM) = IN \\{Dh{z) Ux)- 1 Diag(||z|| d - 1 4 /2 )||, 

or [j,(h, z) = oo if Dh(() | z ± is not invertible. Here, \\h\\ is the Bombieri- 
Weyl norm of h and the second norm in the product is the operator 
norm of that linear operator. Note that fi(h, z) is essentially equal 
to the operator norm of the inverse of the Jacobian Dh((), restricted 
to the orthogonal complement of z. The rest of the factors in this 
definition are normalizing factors which make results look nicer and 
allow projective computations. See [22] for more details. Sometimes \i 
is denoted /i n orm or /i pro j, but we keep the simplest notation here. 

The two following results are versions of Smale's 7-theorem, and 
follow from the study of the condition number in [221 EI]- 

Proposition 1. [6, Prop. 4.1] Let f G Vu) an d let h G Hid) be its 
homogeneous counterpart. Let rj = (771, . . . , r] n ) G C n be a zero of f , 
and let ( = (1, rji, . . . , T) n ) G P(C n+1 ) be the associated zero of h. Let 
x G C n satisfy 

II _ n< 3 ~ yl 

Then, x is an affine approximate zero of f , with associated zero rj. 

Proposition 2. [3J Let ( G P(C n+1 ) be a zero of h G H {d ) and let 
z G P(C n+1 ) be such that 

d R {z, C) < d 3/ 2 J h q > where u o = °- 17586 - 
Then z is an approximate zero of h with associated zero (. 

The following result gives a tool to obtain affine approximate zeros 
from projective ones: 
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Proposition 3. P, Prop. 4.5] Let f G Vu) and let h G Tita) be its 
homogeneous counterpart. Let r\ = (r)i, . . . , r] n ) G C n be a zero of f , 
and let ( = (1,771, . . . ,t] n ) G P(C ra+1 ) be the associated zero of h. Let 
z = (zo, . . . , z n ) G P(C n+1 ) be a projective approximate zero of h with 
associated zero (, such that 

( 3-Vr 

(d R (z,C) < d y2 U ° (ht0 suffices.) 
Let z l = Np(h) l (z), where I G N is such that 

Z>log 2 log 2 (4(l + ||r/|| 2 )). 



dji(z, () < arctan 



Li t s 1 (|,...,|). Then. 



\x — rj\\ < 



d^ii(f,r,Y 

In particular, x l is an affine approximate zero of f with associated zero 
7] by Proposition^ 

Thus, if we have a bound on ||?y|| and a projective approximate zero of 
h with associated zero the projective solution (, we just need to apply 
projective Newton's operator N P (h) a few times |~log 2 log 2 (4(l + 1|?7|| 2 ))] 
to get an affine approximate zero of / with associated zero 77. Here, 
by [A] we mean the smallest integer number greater than A, A G R. 
Thus, a solution to Problem [T] follows from a solution to Problem [2] and 
a control on the norm of the affine solutions of / G V(d)- The latter 
can be done either on per case basis or via a probabilistic argument as 
in P, Cor. 4.9], where it is proved that for / such that ||/|| = 1 and 
5 G (0, 1), we have ||t7|| < Vy/nn/d with probability greater than 1 — 5. 

From now on we center our attention in Problem [2j and we will 
assume that all the input systems h have unit norm, namely h G S. 

3. The homotopy method: A one-root finding algorithm 

Let V = {(/,C) e S x P(C n+1 ) : /(C) = 0} be the so-called solution 
variety. Elements in V are pairs (system, solution). Consider the 
projection on the first coordinate 7r : V — > §. The condition number 
defined above is an upper bound for the norm of the derivative of the 
local inverse of tt near 7r(/, (), see for example [21 Chapter 12]. In 
particular, 7r is locally invertible near (/, Q if //(/, Q < 00. 

Let t — > h t G §, < t < T be a C 1 curve, and let ( be a solution of 
ho- If M^O) Co) < 00, then n is locally invertible near h . Thus, there 
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exists some e > such that for < t < e the zero Co can be continued 
to a zero ( t of h t in such a way that t — > ( t is a C 1 curve. We call the 
curve t — > (h t , ( t ) the lifted curve of t — > h t . There are two possible 
scenarios: 

• Regular: The whole curve t — > h t , < t < T can be lifted to 

(htXt)] 

• Singular: There is some e < T such that t — >■ h t can be lifted 

for < t < e, but [i(h t , Q) — > oo as t — > e. 

Problem 3. Create a homotopy continuation algorithm, a numerical 
procedure that follows closely the lifted curve. Namely, in the regular 
case such algorithm's goal is to construct a sequence — to < ti < . . . < 
tk — T and pairs (gt, Zi) G § x P(C n+1 ) such that for all i = 0, . . . k we 
have gi = h ti and z-i is an approximate zero associated to the zero Q of 
gi with (go, Co) and (gi, Q) lying on the same lifted curve. 

The homotopy method that we have in mind would solve the problem 
above (in the regular case) and would create an infinite sequence {ti} 
converging to the first singularity on the curve in the singular case. 

Remark 1. A homotopy algorithm still may be useful in a singular 
case where the curve can be lifted for t e [0, T), which is the scenario, 
e.g., of a homotopy curve leading to a singular solution. One may use 
Zi for ti close to T as an empirical approximation of the singular zero. 
Approximate zeros (defined before) associated to a singular zero might 
not exist, since Newton's method loses its quadratic convergence near a 
singularity. 

Given a C 1 curve t — > h t , we denote h t = 4 /it- Namely, h t is 
the tangent vector to the curve at t. Note that h t depends on the 
parametrization of the curve, not only on the geometric object (the arc 
defined by the curve). 

A continuous curve t — > h t E S, < t < T is of class C l+Ltp if it is 
of class C 1 in [0,T] (i.e. it has a continuous derivative in (0, T) and 
one-sided derivatives at t — and t = T making h(t) continuous in 
[0, T]), and if the mapping t — > h t is a Lipschitz map, namely if there 
exists a constant K > such that 

\\h t - h s \\ < K\t - s\, Vt,s e [0,T]. 

By Rademacher's Theorem, this implies that the second derivative h t 
exists almost everywhere and is bounded by \\h t \\ < K. 

3.1. Explicit construction of the homotopy method. A certified 
homotopy method and its complexity was shown for the first time in 
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[221 123] . at least for the case of linear homotopy. In a recent work 
|21j . the theoretical complexity of such methods was greatly improved 
although no specific algorithm was shown as the choice of the step 
size was not specified. This last piece can be done in several ways, 
see [3j [lj [10] . We now recall the homotopy method of [3] , designed to 
follow a C l+Lip curve t -> h t G S, t G [0, T\. We assume that: 

(1) We know an approximate zero zo, \\z \\ = 1 of go = h Q , satisfying 

(3.1) d R (z Q ,(a) < ,o /2 U ?, > ? , where u = 0.17586, 

2d- i / 2 fi(h , Co) 

for some exact zero Co of ho. 

(2) Given t G [0, T], we can compute h t and h t — ^r, 

(3) We know some real number H > satisfying 

(3.2) < d 3/2 tf fe|| 2 , 

for almost every t G [0, T]. From now on, we denote 
P = V2 + V4 + 5H 2 G R. 

For z > 1, define (gi+i, z i+ i) inductively as follows. Let a representative 
of Zi be chosen such that ||zj|| = 1. Let s G [0, T] be such that /i s = ^ 
and let & = h s G be the tangent vector to the curve t h t at 
t = s. Let 

(Vdx \ 



(3-3) 



Xi,i 



Dgi(zi 

z* 



r d,, 



(3-4) 



Xi,2 



9i\r + 



V 

2* 



ij 



-1 



o 



1/2 



and consider 

(3-5) 

Let 



<Pi = Xi.\Xi,2- 



;i - V2u /2)^ 



1 



1 + y/2u /2 
and let ti be chosen in such a way that 



1 




(3.6) 

or ti = T — s if 



<U< 



2Pd?/ 2 (p — — s ' ^°t e that this last case happens when 
the step tj chosen with the formula above takes us beyond the limits 



10 



CARLOS BELTRAN AND ANTON LEYKIN 



of the interval [0,T]. The lower bound on (3.6) is used to guarantee 
that the homotopy step is not too small (and thus the total number of 
steps is not too big!). 

Note that in order to compute tpi we must compute the norm of a 
vector (for \i,2) and the norm of a matrix (for Xi,i)- However, we only 
need to do these tasks approximately for we just need to compute a 
number in [ipi,2ipi\. 

In Section 13.31 below we describe the value of the constants to be 
taken in the case of linear homotopy. 

Let g i+ i = h s+u and let 

N ¥ (g i+1 )(zi) 

z i+l 



\N w (g i+1 ){zi)\y 

This way we generate (<7i,zi), (g2,z 2 ), etc. We stop at k such that 
9k = hr, and we output z k e P(C n+1 ). 

3.2. Convergence and complexity of the homotopy method. 

The homotopy method is guaranteed to produce an approximate zero 
of the target system h = if we are in the regular scenario. Moreover, 
its complexity (number of projective Newton's method steps) is also 
well understood and attains the theoretical result of [21J. With the 
notations above, let 

C = [ T ix(h t ,Ct)\\(htXt)\\dt. 



JO 

The reader may observe that Cq is the length of the path (h t ,(t) in the 
condition metric, that is the metric in the solution variety V obtained 
by pointwise multiplying the usual metric (inherited from that of the 
product § x P(C n+1 )) by the condition number \i. 

Theorem 1. [3] With the notations and hypotheses above, assume that 

d R (z , Co) < ,o /2 M °, > \ ? uo = 0.17586. 
2GP/ i /i(/i , Co) 

Then, for every i > 0, Z{ is an approximate zero of gi, with associ- 
ated zero Cij the unique zero of g,i that lies in the lifted path (h t ,(t)- 
Moreover, 

If Co < oo, there exists k > such that hr = 9k- Namely the number 
of homotopy steps is at most k. Moreover, 

k < \Cd 3 / 2 C ] , 
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where 

c= 2P 1 1 l + V2u /2 

In particular, if Cq < oo the algorithm finishes and outputs z^, an 
approximate zero of f — with associated zero (k, the unique zero of 
f that lies in the lifted path (ht,(t)- 

Remark 2. As |~A] < A + 1 for A G IR, we have that the number of 
steps is at most 

1 + Cd 3/2 C . 

Remark 3. // the curve t — » h t is piecewise C 1+Lip we may divide 
the curve in L pieces, each of them of class C 1+Ltp and satisfying a.e. 
\\ht\\ < d 3 / 2 H\\h t || 2 for a suitable H > 0. The algorithm may then be 
applied to each of these pieces and an upper bound on the total number 
of steps is at most 

L + Cd 3 / 2 C . 

Remark 4. If more than one approximate zero of g — fo is known, the 
algorithm described above may be used to follow each of the homotopy 
paths starting at those zeros. From Theorem^ if the approximate zeros 
of g correspond to different exact zeros of g, and if C is finite for all 
the paths (i.e. if the algorithm finishes for every initial input), then 
the exact zeros associated with the output of the algorithm correspond 
to different exact zeros of f = Ht- 

3.3. Linear homotopy. In the case of linear homotopy, the arc-length 
parametrization of the path is 

(3.7) t^h t = g cos(t) + f ~ R l iU ;f\% sin(t), te[0,T], 

l - Re{{f,g)) 2 



where 

T = arcsin \Jl — Re(f, g) 2 = distance^, /) G [0, n]. 

Note that this is a C°° parametrization so in particular it is C l+Lip . 
From (SI Section 2.2], in this case we may take the following value of 
c/P in the description of the algorithm, 

p = 0.04804448... 

The procedure of certified tracking for a linear homotopy is presented 
by Algorithm [TJ 

Algorithm 1. z* = Tr ack Linear Homotopy (f , g, z Q ) 
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Require: f,g G S; Zq is an approximate zero of g satisfying (3.1). 
Ensure: z* is an approximate zero of f associated to the end of the 
homotopy path starting at the zero of g associated to zq and defined 
by the homotopy (3.7). 
1: i «- 0; Si = 0. 
2: while Si ^ T do 
3: 



Compute 

9i <~ h s 



-gsm{s) 



f-Re({f,g))g 
y/1- Re((f,g)y 



: COS S 



I- 

5: 



at s = Si. 

Determine <fi using (3.3), (3.4), and (3.5). 
Let ti be any number satisfying 

0.04804448 0.04804448 

< U < 



6. 

7. 

8 

9 
10. 
11 
12 



if ti > T — s then 

U<-T-s. 
end if 

s i+l ^~ s i + ti; gi + i <— h s . +l ; Zi+x 

i 4- i + 1. 
end while 
z* <— Zt- 



d 3 / 2 tfi 



\\N P (g i+1 )( Zi )\\ 



The bound on the number of steps in Algorithm [T] given by Theo- 
rem [T] is 



(3i 



k < \71d 3/2 C ] . 



4. Finding all roots 

Let us consider polynomial functions in Ou\, where Ou\ is one of 
{'P(d),'H(d),§} with zeros in O n , where O n is either C n or P(C n+1 ). 

Consider a homotopy t — > h t G O^), t G [0,T], connecting the target 
system hx and the start system h along with a set of start solutions 
Z C h \0) C O n . 

Suppose the homotopy curve t —> h t can be lifted to t — > (h t , Ct) G 
{d) x O n , t G [0, T] such that the projection map vr : {d) xOM {d) 
is locally invertible at any t G [0, T). A homotopy path is defined as the 
projection of such lifted curve onto the second coordinate. If the map 
it is locally invertible at t = T as well, then the path is called regular. 

The homotopy t — > ht is called optimal if every £ G Z is the begin- 
ning of a regular homotopy path. If every solution of hr is the (other) 
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end of the homotopy path beginning at some Co £ Zo then we call the 
homotopy total. 

The homotopy method described in Section [3] can also be applied to 
follow all homotopy paths of an optimal homotopy in S. Namely, if 
we have a C l+Lip curve t — )■ h t satisfying the hypotheses of Theorem 
[l] and start with approximate zeros Z associated to a set of solutions 
Zq C Hq 1 (0) C P(C n+1 ), then we may perform the algorithm for each 
of the initial pairs (/to, zq), zq E Zq. By Theorem [TJ the output will be 
approximate zeros associated to §(Zq) distinct solutions of Ht- 

The area of numerical algebraic geometry (see, e.g., [26]) relies on the 
ability to reliably track optimal homotopies and find all roots of a given 
O-dimensional polynomial system of equations in Ou)- To accomplish 
that one has to arrange a total homotopy. 

4.1. Total-degree homotopy. For a target system in E 'P(d), 
(d) = (di, . . . , d n ), define a total degree linear homotopy to be 

(4.1) t^h t = (T-t)ho + ythr, 7 e C*, t E [0,T], 
where the start system is 

(4.2) h = (xf l -l,...,x^-l)eV w . 

There are total degree many, i.e., di---d n , zeros of h that one can 
readily write down. 

Proposition 4. Assume that Ht has a finite number of zeros, and let 
Z be the set of zeros of h above. Then for all but finitely many values 



of the constant 7 the homotopy (4-1) is total. 



If the target system E Vu) has total-degree many solutions, then 
(for a generic 7 ) the homotopy is optimal. 

If the target system E Vu) has fewer than total degree many 
solutions then: 

• some solutions of the target system may be multiple (singular); 

• in case of O n = C n , some of the homotopy paths may diverge 
(to infinity) when approaching t = T. 

Remark 5. To compute singular solutions one may track regular ho- 
motopy paths to t = T — e for a small e > (as in Remark^ and then 
use either singular endgames Section 10.3] or deflation [T71 [TS] . 

To avoid diverging paths one may homogenize the homotopy passing 
from V{d) to H{d) ■ 

Remark 6. Normalizing an optimal homotopy t — )■ h t with respect to 
the Bombieri-Weyl norm brings the system from Hid) to §. In case of a 
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linear homotopy, arc-length re-parametrization enables an application 
of Algorithm^ for each start solution. 

4.2. Other homotopy methods. There are other ways to obtain all 
solutions with homotopy continuation that exploit either sparseness or 
special structure of a given polynomial system, here we list a few: 

• Polyhedral homotopy continuation based on [13] allows to re- 
cover all solutions of a sparse polynomial system in the torus 
(C*) n . 

• In many cases presented with a parametric family of polynomial 
systems it is enough to solve one system given by a generic 
chioice of parameters. Then, given another system in the family, 
the chosen generic system may be used as a start system in 
the so-called coefficient-parameter or cheater's homotopy [26j 
Chapter 7] recovering all solutions of the latter. 

• Special homotopies: e.g., Pieri homotopies coming up in Schu- 
bert calculus [12] are total and optimal by design. 

For the purpose of this paper we assume that some regularization pro- 
cedure (see Remark [5]) has been applied to make these homotopies 
optimal and they are brought to § as in Remark [6] . 

5. Random linear homotopy and polynomial time 

Suppose, given a regular system / e H(d), we would like to construct 
an initial pair (g, ( Q ) in a random fashion so that every root of / is 
equally likely to be at the end of the linear homotopy path determined 
by this initial pair. A simple solution to this problem would be to 



take g to be the start system (4.2) of the total-degree homotopy - 



or its homogenized version - and pick (o from the start solutions with 
uniform probability distribution on the latter. It has been very recently 
proved by Biirgisser and Cucker p] that this is a pretty good candidate 
for the linear homotopy starting pair, as the total average number of 
steps for each path is 0(d 3 Nn d+1 ) that is 0(A^ log(log(Ar)) ), hence close 
to polynomial in the input size, mainly when n » d. 

In [5j El [7J a probabilistic way to choose the initial pair was proposed. 
We now center our attention in the most recent of these works [7J where 
it is proved that, if the initial pair (g,Co) is chosen at random (with 
a certain probability distribution), then the average number of steps 
performed by the algorithm described in Section [3] is 0(d 3 ^ 2 nN), thus 
almost linear in the size of the input. It is also proved that in this way 
we obtain an approximation of a zero of /, so that each of the zeros of 
/ are equiprobable if / has no singular solution. In [H] it is seen that 
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some higher moments (in particular, the variance) of that algorithm 
are also polynomial in the size of the input. In this section we describe 
in detail how the process of randomly choosing (g, (o) works and we 
recall the main results of [7J [8] . 

Given £ G P(C n+1 ), we consider the set 

R c = {heH {d) :h(O = 0,D~h(O = 0}. 

Note that R$ is defined as the set of polynomials in %{d) whose co- 
efficients (in the usual monomial basis) satisfy n 2 + 2n linear homo- 
geneous equalities. Thus, is a vector subspace of Hu). Moreover, 
let e = (1, 0, . . . , 0) T . Then, R eo is the set of polynomial systems 
h = (hi, ... , h n ) G H(d) such that all the coefficients of hi containing 
Xq' or Xq 4-1 are zero, namely 

hi = Xq 1 2 p2,i(Xi, . . . , X n ) + Xq 1 3 p3j(Xi, . . . , X n ) + • • • , 

where pk,i is a homogeneous polynomial of degree k with unknowns 
Xi, . . . , X n . 

Recall that N + 1 is the (complex) dimension of Hm)- The process 
of choosing (g, ( ) at random is as follows: 

(1) Let (M,l) G C n2+n x c N+1 ~ n2 ~ n = C N+1 be chosen at random 
with the uniform distribution in 

B(C N+1 ) = {re C N+1 : \\r\\ 2 < 1}, 

where || • H2 is the usual Euclidean norm in C N+1 . Thus, M is 
a (n 2 + n)-dimensional complex vector, that we consider as a 
n x (n + 1) complex matrix. Note that choosing ||(M, l)\\2 < 1 
implies that ||M||i? < 1 and indeed the expected value of ||M|||i 
is ^rf ■ At this point we can discard / and just keep M. Note 
that this procedure is different from just choosing a random 
matrix, as it induces a certain distribution in the norm of the 
matrix which is precisely the one that we are interested in. 
Hence, choosing (M, I) in the unit ball and then discarding I is 
not a fool job! 

(2) With probability 1, the choice above has produced a matrix 
M whose kernel has complex dimension 1. Let £0 be a unit 
norm element of Ker(M), randomly chosen in Ker(M) with 
the uniform distribution (we may just obtain any such £ by 
solving M( = with our preferred method, and then multiply 
Co by a uniformly chosen random complex number of modulus 
1). Let V be any unitary matrix such that V*( = e . Choose 
a system h at random in the unit ball (for the Bombieri-Weyl 
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norm) of R eo . Then, consider h = h o V* . (This last procedure 
is equivalent to choosing a system at random with the uniform 
distribution in B(R^ ) = {h G R( : \\h\\ < 1}.) 
(3) Let g G Tt(d) be the polynomial system defined by 



g{z) = ^l-\\Mf F h{z) + 
(4) Let 



Mz 



\ {zAv) dn - l Vd~n) 



9 

9 



\9\\ 

Then, we have chosen (g,(o) an d the reader may check that 
9 (Co) — 0, so Co is an exact zero of g. 
Consider the randomized algorithm defined as follows: 

(1) Input / G S 

(2) Choose (g, Co) at random with the process described above 

(3) Consider the path 

t->ht = gcos(t) + — = ==== sm(t), te[0,T\, 

V 1 - Re ((f>9)r 

where T = arcsin a/1 — Re(f, g} 2 , and note that h = g,h T = 
f. Use Algorithm [T] to follow the path h t and output an ap- 
proximate zero of /. 

For given / G §, let NS(f) be the expected number of homotopy steps 
performed by this algorithm, on input / G §. We have seen in (3.8) 
that 

NS(f) < [71rf 3 / 2 C ] 
The main theorems of [S] are now summarized as follows. 

Theorem 2. If f G § is such that every zero of f is non-singular 
(thus, f has exactly V — d\- ■ ■ d n projective zeros), then: 

• The algorithm above finishes with probability 1 on the choice of 
(g,Co), and 

• Every zero of f is equally probable as the exact zero associated 
with the output of the algorithm ( which is an approximate zero 
off)- 

Assuming f G S is chosen at random with the uniform distribution on 
§, the expected value and variance of NS(f) satisfy 

V{NS{f)) < C in Nd 3/2 , Vax(JVS(/)) < C 2 n 2 N 2 d 3 \n(V), 
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where G\ and C2 are universal constants. One may choose C\ = 
71n/V2. 

Note that this theorem not only gives a uniform distribution of the 
probability of producing any given root of a regular system, but also 
gives a good expected running time, with a number of steps which is 
almost linear in the size of the input. 

An algorithm for finding all solutions of a system / with regular 
zeros follows from Theorem |2} repeatedly create and follow random 
homotopies to find one root of the system until total-degree many 
roots are found. Assuming a very small probability of finding less 
than the all the roots, it suffices to choose T> log D such random homo- 
topies. Thus, the expected number of steps of the proposed procedure 
is 0(d?l 2 nNT> logP), which grows fast as the total degree of the sys- 
tem increases. This fast growth is necessary if we are attempting to 
find all the T> solutions of the system. The bound 0(d 3 ^ 2 nNT) logI>) 
is the smallest proven value for the complexity of finding all roots of 
a system. However, this algorithm may not be the most practical one. 



Using the naive start system (4.2) should require, according to p] an 
average number of steps O (d 3 ^ 2 n d+1 NT>) which is a bigger bound that 
0{d?/ 2 nNV logX>), but guarantees that just T> homotopy paths have 
to be followed. 

6. Implementation of the method 

The computer algebra system Macaulay2 - to be more precise, NAG4M2 
(internal name: NumericalAlgebraicGeometry) package [15] - hosts 
the implementation of Algorithm [TJ which is the first implementation 
of certified homotopy tracking in a numerical polynomial homotopy 
continuation software. The current implementation is carried out with 
standard double floating point arithmetic without analyzing effects of 
round-off errors. For a variant of the algorithm that facilitates rigorous 
error control see [I]. 

6.1. NAG4M2: User manual. There are several functions that we 
would like to describe here. First let us give an example of launching 
track procedure with the certified homotopy tracker: 

11 : loadPackage "NumericalAlgebraicGeometry"; 

12 : R = CC[x,y,z] ; 

13 : T = {x~2+y~2-z~2, x*y}; 

14 : (S,xO) = totalDegreeStartSystem T; 

15 : xl = first track(S,T,xO, 

Predictor=>Certif ied , Normalize=>true) 

o5 = {.00000207617, -.706804, .70744} 
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05 : Point 

i6 : xl . NumberOf Steps 

06 = 129 

The values for the optional arguments Predictor and Normalize spec- 
ify that the certified homotopy tracking is performed and the polyno- 
mial systems are normalized to the unit sphere S. In this particular 
example, totalDegreeStartSystem creates an initial pair based on the 
homogenization of the system described in ^2 and track follows the 
linear homotopy starting at this initial pair and finishing at the given 
target system. 

The user can also get a good initial pair (7.1) discussed below with 
the function goodlnitialPair as well as a random pair of start system 
and solution as described in Section [5] with randomlnitialPair. 

It is possible for track to return a solution marked as failure. This 
happens when the step size becomes smaller than the threshold set by 
the optional parameter tStepMin, which has the default value of 10 -6 . 



6.2. Uncertified homotopy continuation. All existing software, such 
as HOM4PS2 [H], Bertini [9J, and PHCpack [27], utilize algorithms 
based on alternating predictor and corrector steps. Here is a summary 
of operations performed at a point of continuation sequence t e [0, T] 
starting with a pair (h t ,x t ) where x t approximates some zero rjt of h t : 

(1) Decide heuristically on the step size At that predictor should 
take; 

(2) Use a predictor method, i.e., one of the methods for numerical 
integration of the system of ODEs 

z = -(Dh t y l h t 

to produce an approximation of (t+At, a solution of h t+ Au 

(3) Apply the corrector: perform a fixed number I of iterations of 
Newton's method to obtain a refined approximation Xt+M = 

N(h t+At ) l (x t+A t); 

(4) If the estimated error bound in step[3]is larger than a predefined 
tolerance, decrease At and go to step [T] 

After tuning the parameters, e.g., tolerances values, the application of 
described heuristics often produces correct solutions. 

6.3. Certified vs Heuristic. We can imagine several "unfortunate" 
scenarios when two distinct homotopy paths come too close to each 
other. Consider sequences z , z tl , . . . , z tk and z' Q , z' t , , . . . , z' t , created by 

an uncertified algorithm in an attempt to approximate these two paths: 
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• If there are subsequences in two sequences that approximate a 
part of the same path then this is referred to as path jumping. 

• Path crossing happens when the sequences jump from one path 
to the other, but there is no common path segment that they 
approximate. 

While path jumping can be detected, in principle, a posteriori and the 
continuation rerun with tighter tolerances and smaller step sizes, the 
path crossing can not be determined easily. 

Path crossing does not result in an incorrect set of target solutions; 
however, for certain homotopy-based algorithms such as numerical ir- 
reducible decomposition [22] and applications relying on monodromy 
computation such as [16] the order of the target solutions is crucial. 
Therefore, one not only needs to certify the end points of homotopy 
paths, but also has to show that the approximating sequences follow 
the same path from start to finish. The certification of the sequence 
produced in Section [3] provided by Theorem [T] gives such guarantee. 

In certain cases the target solutions obtained by means of uncertified 
homotopy continuation can be rigourously certified after all of them are 
obtained. For instance, suppose a target system h T G H{d) has distinct 
regular solutions in P(C n+1 ), then there are total degree many of them. 
Suppose some procedure provides total degree many approximations 
to solutions. If a bound on max{/i(/z T , () \ ( G /i^ 1 (0)} is known, then 
using Proposition [2] these approximations may be certified as distinct 
numerical zeros, thus certifying that all solutions have been found. If 
no such bound is known, one may still try to prove that the zeros are 
different by means of Smale's a-theorem [21] (see [H]). However, these 
procedures do not detect if path crossing has occurred. 

7. Experimental results 

The developed and implemented algorithm provided us with a chance 
to conduct experiments that illuminated several aspects in the complex- 
ity analysis of solving polynomial systems via homotopy continuation. 

7.1. Practicality of certified tracking. Our experiments in this sec- 
tion were designed to explore how practical the certified tracking pro- 
vided by Algorithm [T] is. As was already mentioned, the proposed cer- 
tified procedure makes sense only for a regular homotopy. Moreover, in 
nearly singular examples the certified homotopy is bound to show bad 
performance due to steps being minuscule at the end of paths, which 



is mandated by (3.6). 

In the table below we give the data produced by tracking of total- 
degree homotopy that is optimal for the chosen examples: 
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• Random(rf lj ... j rf n ): a random system in § C H(d) with uniform 
distribution; 

• Katsura n : a classical benchmark with one linear and n — 1 
quadratic equations in n variables. 

For every experiment we provide the number of solutions, the av- 
erage number of steps per homotopy path both for the certified algo- 
rithm (C) and one of the best heuristic procedures (H) implemented in 
Macaulay2. 



system 


#sol. 


#steps/path (C) 


#steps/path (H) 


Random( 2j 2) 


4 


198.5 


31 


RandoHi( 2i 2,2) 


8 


370.125 


23 


Random (2i 2,2,2) 


16 


813.812 


44.375 


RandoHi( 2i 2,2,2,2) 


32 


1542.5 


48.5312 


Random (2i 2,2,2,2,2) 


64 


2211.58 


58.5312 


Katsura 3 


4 


569.5 


25.75 


Katsura4 


8 


1149.88 


41.5 


Katsuras 


16 


1498.38 


39.0625 


Katsura6 


32 


2361.81 


55.5625 



One step in a heuristic algorithm takes more work than that of the 
certified tracker: there is a predictor and several corrector steps per- 
formed and, if unsuccessful, new step size chosen only to repeat the 
procedure. Despite that the heuristic approach leads to much smaller 
computational time for larger systems: it could be concluded from the 
table above that the number of successful heuristic steps does not grow 
fast with degree of the system and the number of variables (in compar- 
ison to certified tracking). 

Of course, it is clear that if we want to run a certified, non heuristic 
method as the one we propose, we will need more computational time. 



7.2. A conjecture by Shub and Smale. In [23], a pair 



(7.1) 







e 



w 



was conjectured to be a good starting pair for the linear homotopy. 
More exactly, let 



Egood = E($(steps) to solve f with lin. homotopy starting at (g,e )), 
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where the expectation is taken for random / G S. Then, the conjecture 
in [23] can be writen a^J 

(7.2) E good < a small quantity, polynomial in N, 

The following experimental data was obtained by running a linear 



homotopy connecting the pair (g, eo) as in (7.1) to a random system 
in § C 'H(d) with di — 2 for i — 1, . . . , n. We compare the values to 
that of B(n,d,N) = IX-nd^^nX j \f2 which according to Theorem [2] is 
a bound for the average number of steps. 



n 


Egood 


Var g ood 


Etotal 


Var t otal 


Brand 


Vciv rand 


B( 


n, 


d,N) 


4 


962.051 


3.2 


10 b 


1263.72 


4.3 


10 b 


1622.29 


6.8 


10 b 


1 





■ 10 5 


5 


1524.6 


6.9 


10 5 


2130.54 


1.2 


10 6 


2728.3 


1.7 


10 6 


2 


3 


• 10 5 


6 


2258.33 


1.3 


10 6 


3129.56 


2.2 


10 6 


4137.16 


3.5 


10 6 


4 


5 


• 10 5 


7 


3130.83 


2.3 


10 6 


4530.55 


4.5 


10 6 


5743.32 


5.5 


10 6 


7 


8 


• 10 5 


8 


4154.38 


3.9 


10 6 


5967.57 


6.7 


10 6 


8048.94 


1.0 


10 7 


1 


2 


• 10 6 


9 


5488.93 


7.0 


10 6 


8013.71 


1.1 


10 7 


10482.1 


1.6 


10 7 


1 


9 


• 10 6 


10 


6871.35 


1.0 


10 7 


10071 


1.4 


10 7 


13477.5 


2.2 


10 7 


2 


9 


• 10 6 


11 


8622 


1.2 


10 7 


12996.1 


2.8 


10 7 


17193.3 


3.5 


10 7 


4 


2 


• 10 6 


12 


10413.3 


2.0 


10 7 


15115.4 


2.8 


10 7 


20761.3 


4.6 


10 7 


5 


8 


• 10 6 


13 


12447.1 


2.6 


10 7 


18744.5 


4.3 


10 7 


25646.5 


6.3 


10 7 


7 


9 


• 10 6 


14 


14769.9 


3.3 


10 7 


22317.1 


6.1 


10 7 


29596.7 


9.1 


10 7 


1 





•10 7 


15 


17255.7 


4.4 


10 7 


26017.7 


7.3 


10 7 


35582.6 


1.2 


10 8 


1 


4 


• 10 7 


16 


20959.7 


5.9 


10 7 


30063.9 


1.0 


10 8 


42098.9 


1.5 


10 8 


1 


7 


•10 7 


17 


23589.4 


7.5 


10 7 


35403.1 


1.3 


10 8 


48024.5 


1.7 


10 8 


2 


2 


•10 7 


18 


27400.9 


9.6 


10 7 


40242.5 


1.5 


10 8 


54955.4 


2.3 


10 8 


2 


7 


•10 7 


19 


29930.3 


1.0 


10 s 


46502.2 


2.3 


10 8 


62855.2 


2.9 


10 8 


3 


4 


• 10 7 


20 


34374.2 


1.4 


10 8 


51730.2 


2.3 


10 8 


71242.5 


3.5 


10 8 


4 


1 


• 10 7 



For each value of n we have generated 1000 random systems in § with 
a uniform probability distribution. The values E goo d and Var goo d are 
estimated expected value and variance of the number of steps taken 
by Algorithm [l] for the initial pair in (7.1); E rand and Var rand refer 
to those for the random initial pair; E tot ai and Var tota i refer to those 
for the homogeneous version of the total-degree homotopy system of 
Section 4.1 containing all the roots of unity (the choice of the root is 
irrelevant for symmetry reasons). Namely, the pair 



(7.3) 



h = (X*-X i 



> 



X. 



in 



x° n 



Co 



The original pair suggested by Shub and Smale had no factors as the one 
here. As done in other papers, we add those factors here to optimize the condition 
number /i(g, eo). 
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The table above and Figure [T] below suggest two conclusions for the 
case of degree two polynomials: 

• The random homotopy seems to take approximately double 



number of steps than the homotopy with initial pair (7.1 ). The 
total degree homotopy lies in-between. 
• The average number of steps in the three cases seem to grow as 
Constant ■ ^= with Constant ~ 35, 50, 70 for E good , E tota i and 
E rand respectively. 

This experiment thus confirms the conjecture by Shub and Smale and 
moreover it suggests a more specific form, suggesting that the same 
bound given for random homotopy should hold for the conjectured 
pair: 

(7.4) E good < CnNd 3 / 2 , 

with C a constant. We also extend this conjecture to the case of the 



initial pair total-degree homotopy pair (h 0} ( ) of Equation (7.3 ) above: 

E to tai < CnNd?' 2 . 

Moreover, as pointed out above, in the case of degree 2 systems, we 
obtained experimentally a better bound 

CN 

Eqoodi Etotai i E ran d ^ . — , 

Jn 



C a constant. The difference between the experimentally observed 
value and the theoretical bound in the case of randomly chosen initial 
pairs, respectively 0(N/y/n) and 0{nN) for (d) = (2, . . . , 2) can be 
explained as follows. The proof of the theoretical bound starts by 
bounding 

Co= f ^{ht,Ct)\\{ht,Ct)\\dt<V2 I fx(h t ,Ct) 2 \\ht\\dt, 



which follows from the fact that \\( t \\ < ^(ht, (t) \\th\\ by the geometric 
interpretation of the condition number. This last inequality is not 
sharp in general, and hence one may expect a better behavior of the 
random linear homotopy method than the one given by the theoretical 
bound. 

7.3. Equiprobable roots via random homotopy. The algorithm 
constructing a random homotopy has been implemented in two vari- 
ants: 

(1) as described in Section [5] 
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(2) the initial pair for the linear homotopy is built by taking (g, eo) 



in ( 7. 1 ) and performing a random unitary coordinate transfor- 
mation (see [19] for a stable and efficient algorithm that chooses 
such a random unitary matrix). 

Then the following experiment was conjured to show the equiprob- 
ability of the roots at the end of a random homotopy promised by 
Theorem [2] as the target system we take / = g + eh where g is as 



in (7.1), h is chosen randomly in S, and e is small. Note that g has 
a unique non-singular solution which is very well-conditioned, but it 
also has a whole subspace of degenerate solutions. Hence, / also has a 
rather well-conditioned solution, and then T> — 1 isolated, but poorly 
conditioned ones. One might expect that the random homotopy (2) we 
have just described (for such a fixed /) would be biased to discover the 
well-conditioned root. Indeed, we obtained numerical evidence that 
this is not the case: all the solutions seem to be equiprobable. 

For / with the degrees d = (2, 2, 2) and e — 0.1 and several random 
choices of g we have made experiments with certified tracking procedure 
making 8000 runs. We experimented with both variants (1) and (2) 
of choosing the random initial pair. Each experiment resulted in close 
to 1000 hits for each of 8 roots — in both variants (1) and (2). This 




Figure 1. In the first figure, we have plotted the ex- 
perimental values obtained for E good , E rand and E tota i 

for n = 4, . . . , 20. In the second one we plot 9 °°^ — , 

g '°'°'" 1/2 and E ^f 2 for n = 4, . . . , 20 
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appears to show the conclusion of Theorem |2j valid for variant (1), and 
moreover extend it to the case of variant (2). 

We can state this experimental result in a more precise way, using 
Shannon's entropy as suggested in [7j. Assume that we have an al- 
gorithm that involves some random choice in its input, and that can 
produce different outputs Xi, . . . ,xi. Shannon's entropy is by definition 
the number 

l 

i=l 

where pi is the probability that the output is x,. It is easy to see that 
Shannon's entropy of an algorithm is maximal, and equal to log 2 (0) if 
and only if every output is equally probable. The experimental value 
of Shannon's entropy for the random algorithm in all experiments de- 
scribed above is in the interval [2.99,3]; the maximum, in this case, is 
log 2 8 = 3. 

The exact reason for the modified algorithm (variant (2)) to produce 
equiprobability of the roots is not understood. This poses a very inter- 



esting mathematical question, which together with proving (7.4) would 
yield a great progress in the understanding of homotopy methods for 
solving systems of polynomial equations. 
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